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NUMERICAL  SOLUTION  OF  THE  NAVIER- STOKES 
EQUATIONS  FOR  2D  HYDROFOILS 


by 

Joe  F.  Thompson,  Samuel  P.  Shanks,  and  Ray  L.  Walker 

Abstract 

This  report  presents  the  results  of  an  investigation  of  the  appli- 
cation of  numerically-generated  boundary-fitted  curvilinear  coordinate 
systems  in  the  finite-difference  solution  of  the  time-dependent,  two- 
dimensional  Navier-Stokes  equations  for  the  laminar  viscous  flow  about 
hydrofoils  moving  either  submerged  at  a finite  depth  or  in  a free  surface 
of  a fluid  of  infinite  depth.  The  hydrofoil  may  be  of  arbitrary  shape, 
and  its  motion  may  include  pitching  oscillation  or  oscillation  normal  or 
parallel  to  the  plane  of  the  undisturbed  free  surface  as  well  as  trans- 
lation parallel  to  this  plane.  A computer  code  has  been  developed  that 
is  capable  of  predicting  the  flow  field,  pressure  distributions,  and 
force  coefficients  for  this  configuration  at  low  Reynolds  numbers.  The 
finite-difference  solution  is  implicit  in  time  so  that  all  the  difference 
equations  are  solved  simultaneously  by  iteration  at  each  time  step. 


I.  INTRODUCTION 


This  report  presents  the  results  of  an  Investigation  of  tin-  appli- 
cation of  numerically- genera ted  houndary-fitted  curvilinear  coordinate 
systems  in  the  finite-difference  solution  of  the  time-dependent,  two- 
dimensional  Navier-Stokes  equations  for  the  laminar  viscous  flow  about 
hydrofoils  moving  either  submerged  at  a finite  depth  or  in  a free  surface 
of  a fluid  of  infinite  depth.  The  hydrofoil  may  be  of  arbitrary  shape, 
and  its  motion  may  include  pitching  oscillation  or  oscillation  normal  or 
parallel  to  the  plane  of  the  undisturbed  free  surface  as  well  as  trans- 
lation parallel  to  this  plane.  A computer  code  has  been  developed  that 
is  capable  of  predicting  the  flow  field,  pressure  distributions,  and 
force  coefficients  for  this  configuration  at  Low  Reynolds  numbers.  The 
finite-difference  solution  is  implicit  in  time  so  that  all  the  difference 
equations  are  solved  simultaneously  by  iteration  at  each  time  step. 

This  investigation  consisted  essentially  of  three  parts: 

(1)  Verification  of  the  use  of  numerically-generated  boundary-fitted 
curvilinear  coordinate  systems  in  Navier-Stokes  solutions  by  application  of 
the  technique  to  t he  unbounded  high  Reynolds  number  flow  past  a semi- inf ini te 
flat  plate  parallel  to  the  undisturbed  flow.  The  results  from  this  solu- 
tion gave  excellent  comparison  with  the  Blasius  boundary  layer  solution 

in  the  region  where  the  latter  is  applicable. 

(2)  Development  of  the  numerical  solution  for  a submerged  hydrofoil 
at  a finite  depth  below  the  free  surface. 

(1)  Development  of  the  numerical  solution  for  a hydrofoil  in  the 
free  surface. 

The  procedures  and  results  of  these  three  solutions  are  summarized  herein. 

The  concomitant  References  1,  2,  and  3,  produced  by  this  investigation. 
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contain  the  details  and  more  extensive  presentations  of  results  from  the 
above-mentioned  Parts  (1),  (2),  and  (3)  of  the  study,  respectively. 

These  three  documents  are  available  from  Mississippi  State  University. 
Results  for  the  finite  flat  plate  of  Part  (1)  have  also  been  reported 
in  References  4 and  5,  while  preliminary  results  from  Parts  (2)  and  (3) 
were  given  in  Reference  6. 

The  numerical  solution  of  the  Navier-Stokes  equations  for  flow  with 
a free  surface  is  complicated  in  particular  by  the  fact  that  part  of  the 
boundary  of  the  calculation  region,  i.e.,  the  free  surface,  is  deforming. 
This  makes  the  accurate  representation  of  boundary  conditions  on  the 
free  surface  difficult;  yet  this  solution,  as  other  partial  differential 
equation  solutions,  is  most  strongly  influenced  by  the  boundary  conditions. 
The  most  critical  need  for  accuracy  thus  lies  in  precisely  the  region  of 
the  most  difficulty  of  attainment. 

Numerical  solutions  for  this  free  surface  flow  problem  have  generally 
tracked  the  moving  free  surface  through  a fixed  grid,  using  interpolation 
among  the  fixed  regularly  spaced  grid  points  to  represent  the  surface 
boundary  conditions.  Similarly,  solid  body  shapes  in  the  flow  have  either 
been  simple,  so  as  to  coincide  with  rectangular  or  cylindrical  grids, 
or  have  been  represented  also  by  interpolation  among  grid  points.  A 
survey  of  such  methods  is  given  in  the  concomitant  Reference  (2). 

The  basis  of  the  present  numerical  solution  is  the  technique  of 
numerically-generated  boundary-fitted  curvilinear  coordinate  systems 
reported  earlier  in  Reference  7.  This  is  a procedure  for  automatic 
numerical  generation  of  curvilinear  coordinate  systems  with  coordinate 
lines  coincident  with  all  boundaries  of  a general  multi-connected. 
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two-dimensional  region  containing  any  number  of  arbitrarily  shaped 
bodies.  The  curvilinear  coordinates  are  generated  as  the  solution  of 
an  elliptic  partial  differential  system.  No  restrictions  are  placed  on 
the  shape  of  the  boundaries,  which  may  even  be  time-dependent,  and  the 
approach  is  not  restricted  in  principle  to  two  dimensions.  With  this 
procedure  the  numerical  solution  of  a partial  differential  system  may  be 
done  on  a fixed  rectangular  field  with  a square  mesh  with  no  interpolation 
required  regardless  of  the  shape  of  the  physical  boundaries,  regardless 
of  the  spacing  of  the  curvilinear  coordinate  lines  in  the  physical  field, 
and  regardless  of  the  movement  of  the  coordinate  system  in  the  physical 
plane.  A number  of  examples  of  coordinate  systems  and  application  thereof 
to  the  solution  of  partial  differential  equations  is  given  in  [8]  and  [9], 
along  with  a discussion  of  the  technique.  This  procedure  essentially  el imi 
nates  the  boundary  geometry  as  a complicating  factor  in  the  numerical  solu- 
tion of  partial  differential  equations.  The  use  of  boundary-fitted  coordi- 
nate systems  for  the  solvit  ion  of  the  incompressible  Navier-Stokes  equations 
for  the  flow  about  two-dimensional  airfoils  has  been  reported  by  Thames, 
et.  al.  [41  and  by  Hodge  [10].  The  latter  reference  uses  the  pressure- 
velocity  formulation  used  in  the  present  work. 

The  use  of  boundary- f i t ted  coordinate  systems  is  particularly 
attractive  for  free  surface  problems,  since  a coordinate  line  will  remain 
coincident  with  the  free  surface  as  it  deforms  under  wave  action.  The 
physical  flow  field  is  transformed  to  the  curvilinear  coordinate  system 
as  discussed  in  more  detail  in  Section  II.  The  field  in  the  transformed 
plane  is  rectangular  with  a fixed  square  grid  regardless  of  the  movement 
of  the  physical  boundaries.  With  the  partial  differential  equations  of 
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motion  and  their  associated  boundary  conditions  transformed  to  the 
curvilinear  system,  all  computation  can  be  done  on  the  fixed  square  grid 
in  the  transformed  plane  regardless  of  the  motion  of  the  free  surface 
or  the  hydrofoil.  It  is  even  possible  to  allow  the  hydrofoil  to 
oscillate  without  really  complicating  the  problem,  since  a coordinate 
line  can  also  remain  coincident  with  the  hydrofoil  surface  as  it  oscillates. 

The  present  solution  is  capable  of  treating  the  low  Reynolds  number 
viscous  flow  about  a translating  hydrofoil  in  or  below  the  free  surface 
of  a fluid  of  infinite  depth.  The  hydrofoil  may  also  be  in  pitching, 
plunging,  or  longitudinal  oscillation  as  well  as  translation.  The  hydro- 
foil starts  from  rest  with  a flat  surface  and  accelerates  to  full  speed 
at  any  acceleration  desired.  The  general  solution  procedure  is  discussed 
in  Sections  II-IV,  and  facets  peculiar  to  each  of  the  three  parts  of  the 
study  together  with  typical  results  are  given  in  Sections  V-VIII. 


4 


II.  BOUNDARY- FITTED  COORDINATE  SYSTEM 


The  basic  idea  of  the  boundary-fitted  coordinate  systems  is  to 
numerically  generate  a curvilinear  coordinate  system  having  some  coordi- 
nate line  coincident  with  each  boundary  of  the  physical  region  of  inter- 
est, regardless  of  the  shape  of  these  boundaries.  This  is  done  by  taking 
the  curvilinear  coordinates  to  be  solutions  of  an  elliptic  partial  dif- 
ferential system, with  constant  values  of  one  of  the  curvilinear  coordi- 
nates specified  as  Dirichlet  boundary  conditions  on  each  boundary.  Values 
of  the  other  coordinate  are  either  specified  in  a monotonic  variation  over 
a boundary  as  Dirichlet  boundary  conditions,  or  are  determined  by  Neumann 
boundary  conditions  thereon.  In  the  latter  case,  the  curvilinear  coordi- 
nate lines  can  be  made  to  intersect  the  boundary  according  to  some  speci- 
fied condition,  such  as  normalcy  or  parallel  to  some  given  direction.  It 
is  also  possible  to  exercise  control  over  the  spacing  of  the  curvilinear 
coordinate  lines  in  the  field  in  order  to  concentrate  lines  in  regions  of 
expected  high  gradients. 

In  any  case,  the  numerical  generation  of  the  coordinate  system  is 
done  automatically  for  any  shape  boundaries,  requiring  only  the  input 
of  points  on  the  boundary.  The  technique  has  been  described  in  detail 
in  earlier  reports  [7-9],  and  the  computer  code,  together  with  instructions 
for  and  examples  of  its  use  in  the  numerical  solution  of  partial 
differential  equations,  is  given  in  [9]. 

The  technique  is  described  in  general  in  this  section.  Each  of 
the  three  parts  of  the  present  study  used  a different  variation  of  the 
basic  procedure  as  is  discussed  for  each  configuration  in  Sections  V-VII. 
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Consider  transforming  the  two-dimensional,  doubly-connected  region  D, 


bounded  by  two,  simple,  closed,  arbitrary  contours,  F^  and  onto  a rectan- 

gular region,  D*,  as  illustrated  in  Fig.  1.  We  require  that  F^  map  onto  F^f, 
r2  onto  r2*,  r3  onto  F^*,  and  onto  . Note  that  and  T2*  are  required 
to  be  constant  q-lines,  while  the  arbitrary  cut  between  contours  F^  and  F2  (i.e., 

r„  and  T.)  becomes  constant  £-lines.  The  region  D is  the  physical  plane,  D* 

3 4 

the  transformed  plane,  the  body  contour,  and  T2  the  remote  boundary  contour. 

As  discussed  previously,  the  curvilinear  coordinates  are  generated  by 
solving  an  elliptic  system  of  suitable  form.  One  such  system  is 

5xx  * ■ P(t>n)  <la) 

nxx+  \y  " Q^’n)  (lb) 

with  Dirichlet  boundary  conditions,  one  coordinate  being  specified  to  be  equal 
to  a constant  on  the  body  and  equal  to  another  constant  on  the  outer  boundary, 
with  the  other  coordinate  varying  monotonically  over  the  same  range  around 

both  the  body  and  the  outer  boundary.  This  system  was  used  for  the  finite 
flat  plate  and  the  submerged  hydrofoil  in  the  present  study. 


Since  it  is  desired  to  perform  all  numerical  computations  in  the  uni- 
form rectangular  transformed  plane,  the  dependent  and  independent  variables 


must  be  interchanged  in  Eq . (1).  This  results  in  the  coupled  syster 


ax^  - 2Bx^  + YXnr|  “ - J Ix^PU.n)  <-  x^QU.n)] 
ay^  - 26y?Ti  + yynn  = - J2[y^P(^»n)  + ynQ(S,n)] 


where 


2 2 

x + y 
n n 


2.2 

T -t  + 


8 ’ Vn  + Vn 


J ' Vn  * Vt 


i 

! 

I 
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The  system  described  by  Eq.  (2)  Is  a quasl-llnear  elliptic  system  for 
the  coordinate  functions  x(£,n)  and  y(£,n)  in  the  transformed  plane.  This 
set  is  considerably  more  complex  than  the  linear  system  specified  by  Eq.  (1), 
but  the  boundary  conditions  are  specified  on  straight  boundaries,  and  the  coor- 
dinate spacing  in  the  transformed  plane  is  uniform. 


The  £=constant  lines  may  be  spaced  as  desired  around  the  boundaries, 
since  the  assignment  of  the  ^-values  to  the  [x,y]  boundary  points  is 
arbitrary.  (Numerically,  the  discrete  boundary  values  [ , y ^ ] are  trans- 
formed to  equi-spaced  discrete  (j^-points  on  both  boundaries.)  Control  of 
the  radial  spacing  of  the  p=constant  lines  and  of  the  incidence  angle  of 
the  £=constant  lines  at  the  boundaries  is  accomplished  by  varying  the 
functions  P(£,n)  and  Q(£,n)  in  (2).  All  numerical  computations,  both  to 
generate  the  boundary-f xtted  coordinate  system  and  subsequently  to  utilize 
the  coordinates  for  solving  a set  of  partial  differential  equations,  are 
executed  on  a rectangular  field  with  a uniform  grid. 

The  effect  of  changing  the  functions  P(£,n)  and  Q(£,n)  on  the 
coordinate  system  is  discussed  in  Ref.  9.  One  particularly  effective  pro- 
cedure, used  here  for  the  finite  flat  plate  and  submerged  hydrofoil  solu- 
tions, is  to  choose  P and  Q as  exponential  terms,  so  that  the  coordinates 
are  generated  as  the  solutions  of 


5 + S 

xx  yy 


n 

-l  a sgn(C  - E1)exp(-c  |C 


- I b,  sgn(£  - £j)exp(-dj/(£  - + (n  - h^)*) 


£ P(£,n) 


(3a) 
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nxx  + nyy  " " z ai  88n(ri  - n1)e*p(-ci|n  - n± I ) 


m 


E b sgn(n  - nJexp(-d  /(C  - S.)^  + (n  - n.)^) 
J J j j j 


= QU.n) 


(3b) 


where  the  positive  amplitudes  and  decay  factors  are  not  necessarily  the  same 
in  the  two  equations.  Here  the  first  terms  have  the  effect  of  attracting 
the  E,  = constant  lines  to  the  £ = lines  in  Equation  (3a),  and  attracting 
n = constant  lines  to  the  n = lines  in  Equation  (3b).  The  second  terms 
cause  E,  = constant  lines  to  be  attracted  to  the  points  (£j,rij)  in  (3a), 
with  similar  effect  on  q = constant  lines  in  (3b).  Several  examples  of  the 
use  of  coordinate  system  control  are  given  in  Ref.  9. 

The  transformation  technique  developed  above  can  be  used  to  solve 
time-dependent  problems  with  moving  boundaries  by  performing  all  numer- 
ical calculations,  without  interpolation,  on  a fixed  rectangular  field 
with  a uniform  square  grid  in  the  transformed  plane. 

As  discussed  previously,  the  physical  plane  grid  system  is  gener- 
ated by  solving  the  set  of  elliptic  partial  differential  equations, 

(2),  with  one  of  the  (£,q)  coordinates  specified  to  be  constant  on  the 
boundaries  of  the  physical  plane,  and  the  other  (£,q)  coordinate  dis- 
tributed along  the  boundaries  as  desired.  If  the  boundary  values  of  x 
and  y are  changed  in  the  physical  plane  by  the  movement  of  the  free  sur- 
face contours  a new  solution  of  the  elliptic  system  with  the  changed 
boundary  values  is  obtained  over  the  same  range  of  values  of  £ and  n in 
the  field.  Thus,  the  transformed  plane  remains  unchanged  as  the  coordinate 
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grid  system  moves  in  the  physical  plane.  Only  the  values  of  the  physical 
coordinates  (x,y)  change  with  time  at  the  fixed  grid  points  in  the  trans- 
formed plane. 


The  transformed  time  derivative  is 


li-v  = a(*.y»f)  / 3(x..y»t)  = ,11% 

3t'  au.n.t)  3(C,n,t)  V3t' 

x > y s»n 


1_  ._3f  jjjr  _3_f  3y.  . 3x. 

J C3t  3n  “ 3n  3£H3t ’ r 

£»n 


A /-Af.  .lx  H lx-,  / 3y~, 

J C3£  3n  " 3n  3£H3t'_ 

C.n 


(4) 


All  derivatives  are  expressed  in  the  transformed  variables  (£,n); 
thus  eliminating  the  need  for  interpolation  between  points  in  the 
physical  plane.  The  movement  of  the  physical  plane  grid  points  is 
accounted  for  by  the  time  rate  of  change  of  x and  y,  (-ir~)  and 

O L 

C»n 

3y 

in  the  above  expression. 


III.  EQUATIONS  OF  MOTION 


The  equations  of  motion  are  the  complete  time-dependent  Navier- 
Stokes  equations  with  the  gravity  term  included.  The  no-slip  boundary 
condition  is  applied  on  the  hydrofoil,  and  the  viscous  stress  condi- 
tions are  applied  on  the  free  surface.  The  free  surface  deforms  in  time 
as  waves  are  formed  thereon. 

All  quantities  are  non-dimensionalized  with  respect  to  the  trans- 
lation velocity  of  the  hydrofoil  and  the  hydrofoil  chord.  The  Reynolds 
and  Froude  numbers  are  defined  in  terms  of  these  reference  values.  The 
physical  coordinate  system  is  taken  to  be  fixed  relative  to  the  trans- 
lational motion  of  the  hydrofoil.  In  the  physical  plane  the  equations 
of  motion  are 


ut  + (u2)x 

+ (uv) 

y 

= -p  + (u  + u )/R-ty 

rx  xx  yy 

(5a) 

Vt  + (uv)x 

+ (v2)y 

= - p +(v  +v  )/R-  1/F2 

ry  xx  yy 

(5b) 

P + D * - u2  - 

2u  v - v2  - D 

(5c) 

*xx  * yy 

X 

y x y t 

Voc 

V 

where  R = and  F = 

o 

are 

the  Reynolds  and  Froude  numbers, 

respec- 

V 

/gc 

tively,  Vq  being  the  magnitude  of  a reference  hydrofoil  translational 
velocity,  V the  instantaneous  velocity,  c the  chord,  v the  kinematic 
viscosity,  and  g the  acceleration  of  gravity.  The  third  of  these 
equations  is  the  Poisson  equation  for  the  pressure,  derived  by  taking 
the  divergence  of  the  Navier-Stokes  equations  and  requiring  that  the 
continuity  equation  (D  = ux  + v^  = 0)  be  satisfied.  The  time  derivative 
of  D,  ideally  zero,  has  been  retained  in  this  equation  as  a corrective 
term  in  the  manner  of  Hirt  and  Harlow  [11]. 
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The  boundary  conditions  are  as  follows: 

(a)  On  the  hydrofoil  (no-slip  condition): 


u *=  uR(x,y,t) 


v = Vg(x,y, t) 


where  (uD,vD)  are  the  velocity  components  of  the  hydrofoil  surface  at 

O D 

(x,y,t)  relative  to  the  coordinate  system  translating  with  the  hydro- 
foil. (These  values  are  zero  if  the  hydrofoil  is  not  oscillating.) 

(b)  On  the  surface  (viscous  stress  condition): 


I Vl  + R (uy  + Vx)n2  = (P  - Po)nl 


I vyn2  + i (uy  + Vx)nl  = " p0)n2  (6b) 

where  pQ  is  the  applied  pressure  from  the  atmosphere,  and  n^  and  n2  are 
the  components  of  the  unit  normal  to  the  surface  in  the  x and  y direc- 
tions, respectively.  These  relations  assume  no  wind  stress  on  the  sur- 
face. 

(c)  On  the  remote  boundary  (undisturbed  flow) : 


u = - V 


v ■ 0 


p " po  + (yo  “ y)/F 


(These  conditions  apply  on  the  remote  boundary  strictly  only  until  sur- 
face waves  reach  it.  At  low  R,  sufficient  viscous  dissipation  is  pre- 


sent to  damp  the  waves  before  the  remote  boundary  10  chords  distant  is 


reached. ) 


Using  the  derivative  transformation  relations  given  in  the  appen- 
dices of  either  [2]  or  [9],  these  equations  may  be  transformed  to  the 
curvilinear  coordinate  system,  so  that  the  equations  of  motion  in  the 
transformed  plane  are 


ut  - xt(Vs  " ycun)/J  ' yt(Vn  " xnV/J 


+ [yn(u2)^  - y^(u2)fi]/j  + [x^(uv)^  - (uv) ^ ]/J 


+ (ynpC  " yCpn)/J  = (au£5  " 26%  + yunn 


+ ou  + Tu  )/RJ2  + V 

n K 


(7a) 


vt  “ xt(ynvC  " yCvn)/J  - yt(Vn  ' xnV/J 

+ [yn(uv)c  - y^CuvJ^J/J  + [x^(v2)^  - xri(v2)^]/J 

+ (xcpn  ' Xnp£)/J  = (avC£  ■ 26vCn  + yvnn 
+ ov^  + tv^)/RJ2  - 1/F2  (7b) 

ap«  ‘ 2fipCn  + Tpnn  + <’pn  + Tpf,  ' ' (V{  ‘ Vn)2 

- 2<xc\  ' VE>  (V?  ‘ P{V 

- (xC\  ‘ xnV2  - j2|>t 

+ Wn  - V{)J  - Wn  - DnVJ  (7c> 

where 

D = (ynuC  " yCun  + Vn  “ xnV/J  (7d) 
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and 


a i J*Q(£,n)  , T = J2P(£,n)  . 

The  coefficients  a,  0,  y»  and  J have  been  defined  in  the  previous  sec- 
tion. 

The  time  derivatives  have  also  been  transformed  in  these  equations. 
Thus,  time  derivatives  in  Eq.  (7)  are  taken  with  £ and  n fixed,  while 
those  in  Eq.  (5)  were  taken  with  x and  y fixed.  This  transformation  of 
time  derivatives  allows  the  computation  to  be  done  on  a fixed  grid  in 
the  transformed  plane  even  though  the  physical  grid  is  in  motion  due  to 
the  free  surface  and  hydrofoil  movement.  The  same  procedure  was  used  by 
Shanks  f 2]  and  is  discussed  in  more  detail  therein,  as  well  as  in  [9 ] . 
The  transformed  boundary  conditions  are 

(a)  On  the  hydrofoil  (no-slip  conditions): 

u = 

v = VB (£ , t ) 

(b)  on  the  free  surface  (viscous  stress  conditions) : 

uc  = ^ [(“e  ' JVn)un  + (Jxn)vn  + “T  Vp  " po>  1 (8a) 

v£  = ^ [(-  Jy^)un  + (“6  + JVn)vn  " HT  xn(p  " po)]  (8b) 

(c)  On  the  remote  boundary  (undisturbed  flow) : 

u = - V 
v = 0 

p “ pg  + (yQ  - y)/F? 


13 


The  two  free  surface  boundary  conditions  given  above  result  from 

transforming  the  two  viscous  stress  equations  (6)  to  the  curvilinear 

coordinate  system  and  solving  these  two  simultaneous  equations  for  Uj. 

and  v,  in  terms  of  u and  v , the  free  surface  being  a line  of  constant 
5 n n 

£ in  the  configuration  used.  (See  the  appendix  of  [2]  for  this 
development . ) 

On  the  hydrofoil  contour  and  the  free  surface,  the  pressure  is 
determined  by  iteratively  adjusting  the  pressure  at  each  point  on  the 
hydrofoil  in  proportion  to  the  divergence  of  the  velocity  at  the  same 
point,  so  that  upon  convergence  the  continuity  equation  is  satisfied  at 
the  hydrofoil  surface.  Thus  on  the  hydrofoil  surface,  since  u^  = v^  = 0 
by  the  no-slip  condition,  we  have,  using  (7d), 

p(k+l)  = p (k)  _ k (x  v - y u )/J  (9a) 

s n s n 

while  on  the  free  surface. 


p(k+D  = p(k) 


KD 


(9b) 


with  D given  by  (7d).  Here  (k)  is  the  iteration  counter,  and  K is  a 
proportionality  factor  given  by 

K = 

(a  + y)At 


on  the  hydrofoil  and  by 


on  the  free  surface,  u>  being  an  acceleration  parameter.  The  different 
form  on  the  free  surface  results  from  the  need  to  prevent  positive 
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feedback  from  the  surface  stress  condition  to  the  surface  pressure 
iteration.  See  Appendix  D of  [3]  for  the  development  of  these  relations. 


The  y coordinate  on  the  free  surface  is  determined  at  each  time 
from  the  movement  of  the  free  surface.  Since  the  free  surface  can  be 
described  by  y = f(x,t)  or  f(x,t)  - y = 0,  the  convective  derivative 
of  the  latter  function  must  vanish: 


IV.  NUMERICAL  SOLUTION 


All  space  derivatives  in  the  field  are  approximated  by  second- 
order,  central  difference  expressions.  (A5  and  An  are  both  unity  by 
construction,  the  actual  values  of  5 and  n being  immaterial  since 
cancellation  occurs  after  substitution  in  the  transformed  equations.): 


(Vij  = 2 (fi+l,j  ' fi-l,j} 


(Vij  = 2 (fi,j+l  " fi,j-l) 


(frr)J4  = f,.,  . - 2f . . + f,  . . 
55  ij  i+l,J  ij  i-l.J 


(f  ^ f . . - 2f. . + f,  . . 

nn  ij  1,3+1  ij  i, j-1 


(f5n)ij  ' A (fi+l,j+l  " fi+l, j-1  " fi-l, j+1  + fi-l,J-l) 


Derivatives  along  coordinate  lines  emanating  from  the  hydrofoil 
surface  or  from  the  free  surface  are  evaluated  using  second-order,  one- 
sided difference  expressions  of  the  form 


(Vl.j  ■ 2 + 4f2,j  - 3£ia> 


(Vi,l  = 2 ( f i,  3 + Afi,2  " 3fi,l) 


Finally  all  the  time  derivatives  are  approximated  by  first-order, 
backward  difference  expressions,  so  that  the  solution  is  implicit  in 
time.  The  set  of  five  simultaneous  difference  equations  from  (2)  and 
(7),  three  equations  of  motion  and  two  coordinate  system  equations,  with 
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the  boundary  conditions  are  solved  at  each  time  step  by  SOR  iteration. 
The  result  from  the  previous  time  step  serves  as  the  initial  guess  for 
the  iteration  at  the  next.  The  solution  starts  from  rest  with  a flat 
free  surface  and  proceeds  with  a specified  acceleration  to  full  speed. 

The  body  force  components  are  obtained  from  the  integration  of  the 
pressure  and  shear  forces  around  the  wetted  portion  of  the  hydrofoil 
surface: 

Fx  = ~ 2 $ + | $ (11a) 

Fy  = + 2 f PX^  + | 0 (lib) 


with  vorticitv,  u,  given  by 

to  = ( y v„  - y.v  - x.u  + x u_)/J 

;n  ( C n £ n n £ 


(12) 


Here  the  n-derivatives  are  evaluated  by  the  second-order,  one-sided  dif- 
ference expressions  given  above,  while  the  second-order  central  expres- 
sions are  used  for  the  ^-derivatives . 

Finally,  the  lift  and  drag  coefficients  are  given  by 

Cr  «=  F cose  - F sine  (13a) 

L y x 

C_  = F sine  + F cosO  (13b) 

D y x 

where  6 is  the  angle  of  attack. 
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V.  SEMI-INFINITE  FLAT  PLATE  SOLUTION  (Ref.  [1]) 


This  solution  was  developed  in  the  stream  function-vorticity  formu- 
lation of  the  Navier-Stokes  equations  following  Thames  [12],  and  has  also 
been  reported  in  [5]. 

The  stream  function-vorticity  formulation  of  the  two-dimensional, 
incompressible  viscous  flow  equations  is  given  by 


w + A 0)  - «(o)  +o)  )/R 

t y x x y xx  yy 


(14a) 


ill  + ill  - -0) 

xx  yy 


(14b) 


where  is  the  nop-dimensional  stream  function,  to  the  non-dimensional 
vorticity,  and  R is  the  Reynolds  number  based  on  the  characteristic  velocity 
(free  stream  value)  and  body  length.  The  set  (14)  is  in  the  non-conservative 
formulation.  Eqs.  (14)  may  be  transformed  utilizing  the  operations  given 
in  [9]  yielding  the  set  applicable  in  the  rectangular  transformed  plane. 

The  transformed  equations  are: 


u»t  + (^  - i^a>n)/J  - (au.^  - 26^  + yu^/j2*.  + (Qoo^  + Pc^)/R 
(aiji^  - 26^^  + Yi|^)/J^  + QiJ^  + Pij^  = -to 


(15a) 

(15b) 


with  boundary  conditions 


* “ " constant,  /y~ *n/J  - 0,  e r * (15c) 

* - y(4,n2)cos0  - x(4,n2)sin0,  to  - 0,  [e,n2l  e r2*  d) 
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where  0 is  the  angle  of  attack.  The  second  of  (15c)  guarantees  that  the  no- 
slip condition  is  satisfied  on  the  body  surface.  The  satisfaction  of  this 
condition  is  accomplished  by  iteratively  adjusting  the  value  of  the  vorticity 
on  the  body  surface,  utilizing  a false-position  procedure,  until  the  second- 
order  forward  difference  approximation  to  the  velocity  component  tangential  to 

the  body  surface,  VT  = i(/  /J,  is  below  some  tolerance.  The  iterative  algorithm 

~h 

is  given  by 


(k+1 ) 

1,1 


o(k) 

1,1 


.00 


(k— 1 ) 


ui,i  ~ Mi,l 


(k-1) 


(V  )^  - (V  ) 

' I 'i,l  ^ t 'i,l 
*n  ~n 


(v  )^ 

V1*1 


(16) 


where  k denotes  iteration  count,  5 an  adjustable  parameter,  and  (i,l)  refers  to 
some  point  on  the  body  surface.  This  method  is  an  extension  of  an  approach 
suggested  by  Israeli  [13). 

The  transformation  from  the  physical  to  the  transformed  field  is 
indicated  schematically  in  Fig.  2.  The  coordinate  system  was  generated 
as  the  solution  of  (2)  with  boundary  conditions  as  follows: 

on  a'b'  (plate  surface):  x = specified  as  desired,  y = 0. 

h 

on  a'c'  (upper  boundary):  x = specified  as  desired,  y = 10(£)  . 

R 

on  b'c'  (downstream  boundary):  x = constant,  y = specified  as  desired. 

on  a 'a'  (leading  edge):  x = y = 0. 
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The  condition  on  y on  the  upper  boundary,  a'c',  places  this  boundary  at  twice 
the  Blasius  boundary  layer  thickness  [14].  The  downstream  boundary  was  located 
at  multiples  of  the  distance  at  which  the  slope  of  the  Blasuis  boundary  layer 
is  0.01. 

The  boundary  conditions  for  the  Navier-Stokes  equations  in  the  transformed 
plane  are  as  follows: 

on  a'b*  (plate  surface):  = 0 (v  = u = 0,  no-slip  condition), 

on  a'c'  (upper  boundary);  (J  + , u>  = 0(u=l,w  = 0, 

free  stream  conditions). 


on  b'c'  (downstream  boundary): 


U) 


erf  n dn, 

Rv^.  ,.  . . 

-^-)  (infinite  plate 

solution,  Schlichting 
[14]) 


on  a'a'  (leading  edge):  ip  = u>  = 0. 


The  condition  on  the  downstream  boundary,  b'c',  is  the  exact  solution  of  the 

Navier-Stokes  equations  for  a suddenly  accelerated  fully  infinite  flat  plate. 

The  numerical  quadrature  was  done  by  trapezoidal  integration.  The  condition  on 

on  the  upper  boundary  expresses  u = 1,  the  free  stream  velocity.  All  these 

boundary  conditions  were  implemented  directly  except  the  p = 0 condition  on  the 

n 

plate,  a'b',  which  was  satisfied  by  adjusting  the  value  of  the  vorticity  at  each 
point  on  the  body  by  the  false-pos j tion  iteration  procedure  discussed  above. 

The  coordinate  system  used  for  the  semi-infinite  plate,  shown  in  Tig. 1, 
has  a curved  boundary  located  at  twice  the  Blasuis  boundary  layer  thickness 
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above  the  plate,  with  coordinate  lines  coming  to  a point  at  the  leading  edge. 
This  form  was  chosen  in  preference  to  systems  with  rectangular  boundaries  in 
order  to  concentrate  the  coordinate  lines  near  the  plate  to  a greater  degree 
as  the  Reynolds  number  increases  and  also  to  ensure  a test  of  a representa- 
tive non-orthogonal  curvilinear  system. 

Velocity  profiles  obtained  using  this  coordinate  system  are  shown  in 
Fig.  4 and  compared  therein  with  the  Blasius  boundary  layer  solution 
(Schlichting  [14]).  (Positions  are  given  in  fractions  of  the  distance  to  the 
downstream  boundary.)  Since  the  downstream  boundary  condition  was  the  time- 
dependent  solution  for  the  completely  infinite  plate,  for  which  the  boundary 
layer  thickness  increases  without  bound  as  time  increases,  the  agreement  with 
the  Blasius  solution  deteriorates  as  expected  as  this  boundary  is  approached 
(position  1.0  in  these  figures).  The  loss  of  flow  in  the  lower  portion  of 
the  boundary  layer  that  results  from  this  continual  thickening  of  the 
boundary  layer  on  the  downstream  boundary  causes  the  over-shoot  of  the  Blasius 
profile  that  occurs  upstream  of  this  boundary.  The  agreement  with  the 
Blasius  profile  in  regions  farther  removed  from  the  downstream  boundary  is 
good.  Note  that  the  profiles  upstream  cling  to  the  Blasius  as  the  down- 
stream profile  moves  away. 

Coordinate  system  control  was  used  to  cause  the  system  to  expand  down 
the  plate.  The  agreement  with  the  Blasius  solution  extends  very  near  the 
leading  edge,  since  the  coordinate  lines  are  more  closely  spaced  near  the 
leading  edge.  With  the  Blasius  boundary  layer  solution  as  the  downstream 
boundary  condition  the  problem  becomes  a steady-state  problem,  and  the 
agreement  with  the  boundary  layer  solution  is  excellent ,/  except  at  the  first 
few  steps  near  the  leading  edge  [1],  [5], 


VI.  SUBMERGED  HYDROFOIL  SOLUTION  (Ref.  [2]) 


Figure  5 shows  the  basic  doubly-connected  transformation  with  a free 
surface.  This  type  of  configuration  has  been  used  successfully  for  air- 
foils in  previous  studies  [4].  For  a free  surface  problem  would  be  the 
arbitrary  hydrofoil,  C ^ would  be  the  "infinity"  boundary,  and  would  be 
the  free  surface.  Since  the  "infinity"  boundary  is  chosen  to  be  ten  chords 
from  the  hydrofoil  in  the  present  research,  the  contour  C,.  would  be  approxi- 
mately twenty  chords  long.  Thus,  fewer  points  would  be  on  to  cover  20 
chord  lengths  than  would  be  on  to  cover  approximately  2 chord  lengths. 
Unless  many  E;-points  were  used  the  wide  grid  spacing  on  the  free  surface 
would  cause  large  truncation  error. 

Several  modified  coordinate  systems  were  investigated  in  order  to 
provide  more  points  on  the  free  surface  [2]. 

Figure  6 shows  the  transformation  that  was  used  in  this  research. 

This  coordinate  transformation  was  chosen  because  the  number  of  free 
surface  (C^  and  Cg)  grid  points  is  independent  of  the  number  of  points 
on  the  body  (C^) . Also,  no  points  with  zero  Jacobians  occur  on  the  free 
surface . 

The  transformed  plane  of  Figure  6 forms  a T-shaped  region.  The 
lower  part  of  the  coordinate  system  is  the  same  as  the  basic  transforma- 
tion of  Figure  5.  The  cut  C*  is  taken  at  II  - (12  + -i)  , thus  creating 

the  two  common  reentrant  boundaries  C*  and  C*.  The  upper  part  of  the 
transformed  plane  bounded  by  the  constant  n-line  (JM  - y)  C , C^, 
is  added  to  the  basic  transformed  plane  to  provide  free  surface  and 
"infinity"  boundaries.  Ttie  two  common  reentrant  boundaries  and  C,- 
are  created  to  provide  more  points  on  the  infinity  boundary  than  are  on 
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the  body  (C£) . The  cuts  are  taken  at  half  indices  because  the  point 
(II  - JM  - -|)  has  a zero  Jacobian.  By  taking  the  cut  at  half 
indices,  the  zero  Jacobian  point  is  eliminated  from  the  field  calcula- 
tions . 

The  system  of  finite  difference  equations  is  solved  simultaneously 
by  the  successive-over- relaxation  (SOR)  iterative  method.  The  number 
of  simultaneous  equations  to  be  solved  is  (JMAX  - JM  + 1)(IMAX  - 1)  + 

(JM  - 1) (12  - II  + 1).  Boundary  values  are  specified  on  J = JMAX  for 

all  ie[l,  IMAX],  Also,  boundary  values  are  specified  on  j - JMAX  for 

all  ie[Il,  12].  Boundary  values  or  the  Neumann  boundary  condition  x ^ = 

0 (normal  q-lines  to  free  surface)  may  be  expressed  on  the  free  surface 
contours,  i = 1 and  i = IMAX.  At  the  branch  cut  for  the  constant  n-line 
J = JM  and  it.  (1,  II  - 1],  we  have 

(i,  j - 1)  = (IMAX  - i + 1,  JM)  . 

Also,  at  the  branch  cut  for  the  constant  n-line  J = JM  and  ie[I2  + 1, 
IMAX],  we  have 


(i,  j - 1 ) = (IMAX  - i + 1 , JM)  . 

At  the  branch  cut  for  the  constant  £-line  i = II  and  j e [ 1 , JM  - 1],  we 
have 


i - 1 = 12  . 

At  the  branch  cut  for  the  constant  £-]ine  i = 12  and  jc[l,  JM  - 1],  we 
have 


i + 1 = II  . 
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The  basic  hydrofoil  geometry  and  coordinates  are  shown  in  Figure  7. 


L is  the  chord  length  of  the  hydrofoil  used,  h is  the  depth  of  the 
water  from  the  bottom  to  the  undisturbed  free  surface,  d is  the  depth 
of  the  hydrofoil  below  the  free  surface.  U is  the  reference  velocitv 
for  the  problem  (usually  the  steady-state  free  stream  velocitv),  v is 
the  kinematic  viscosity,  p is  the  density,  D is  the  incompressible 
continuity  equation,  and  g is  the  local  gravitational  constant. 

The  equations  of  motion  are  those  given  in  Section  III  except 
that  the  non-conservative  form  of  the  convective  terms  in  the  Navier- 
Stokes  equations  was  used  for  this  solution. 

Typical  results  of  the  numerical  solution  are  given  in  Figures 
8-18  for  a Karman-Tref f tz  hydrofoil  (Figure  8)  and  in  Figures  19-20 
for  a circular  cylinder  hydrofoil.  The  airfoil  is  defined  by  37 
coordinate  points  and  is  located  one  chord  below  the  free  surface. 

The  field  sizes  of  the  coordinate  grid  are  54  x 30  and  54  x 60.  The 
54  x 60  field  has  its  "outer"  boundary  located  20  chords  from  the  air- 
foil, and  the  54  x 30  field  has  its  "outer"  boundary  located  10  chords 
from  the  airfoil.  Using  Figure  6,  JM  =6,  II  = 10  and  12  = 45.  Six  n- 
lines  were  attracted  to  the  airfoil  with  an  amplitude  of  1000  and  a 
decay  factor  of  1.0. 

Figure  9 shows  the  effect  of  Froude  number  on  the  free  surface 
movement.  Three  Froude  numbers  of  0.5,  1.0  and  2.0  are  shown  for  a 
constant  Reynolds  number  of  20  and  at  a time  of  8.0.  All  three  cases 
were  accelerated  gradually  over  a time  of  4.0.  Also,  the  same  time 
step  size  (At  = .01)  was  used  in  all  three  cases.  At  a time  of  8.0  the 
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airfoil  had  moved  4 chords  at  the  free  stream  velocity  of  one.  By  com- 


paring Figures  9a,  b,  and  c,  the  effect  of  the  airfoil  on  the  free  sur- 
face increases  as  the  Froude  number  increases.  This  result  is  to  be 
expected  since  the  Froude  number  is  the  ratio  of  inertial  forces  to 
gravitational  forces.  Therefore,  as  the  Froude  number  increases  the 
inertial  forces  increase. 

Figures  10  , 11  , and  12  demonstrate  the  effect  of  Froude  number  on 
the  drag,  lift,  and  pressure.  The  time  histories  (Figure  10)  of  drag 
for  the  three  Froude  numbers  of  0.5,  1.0  and  2.0  are  presented  for  Re  = 
20  and  t = 8.0.  Referring  to  the  peak  drag  for  each  figure,  the  drag  is 
reduced  as  the  Froude  number  is  increased,  because  as  the  Froude  number 
increases  the  free  surface  rises  over  the  airfoil  which  changes  the 
local  angle  of  attack.  From  Figure  11,  the  lift  changes  drastically 
because  as  the  Froude  number  decreases  the  buoyancy  forces  become  more 
dominant.  From  Figure  12  the  effect  of  Froude  number  on  the  pressure 
distribution  can  be  seen.  Buoyancy  forces  are  dominant  at  a Froude 
number  of  0.5  (Figure  12b),  and  inertial  forces  are  dominant  at  a Froude 

number  of  2.0  (Figure  12d) . 

Figure  13  shows  three  time  slices  of  a coordinate  system.  The  air- 
foil is  accelerated  to  the  free  stream  velocity  over  a time  period  of  4. 
The  flow  parameters  are  F = 1 .0  and  Re  = 20 , and  the  airfoil  is  located 
one  chord  below  the  free  surface.  At  t = 2.0,  the  airfoil  has  achieved 
half  its  free  stream  velocity.  The  airfoil  movement  has  caused  the  free 
surface  to  rise  slightly  in  front  of  the  airfoil,  and  to  sink  slightly 
behind  the  airfoil.  At  t = 4.0,  the  airfoil  has  achieved  the  free 
stream  velocity  of  1.0.  The  free  surface  peak  has  moved  to  about  the 
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quarter  chord  of  the  airfoil,  and  the  lowest  point  on  the  free  surface 
has  moved  one  chord  behind  the  airfoil.  At  t = 6.0,  the  airfoil  has 
moved  2 chords  at  a constant  velocity  of  1.0.  The  first  peak  has 
remained  at  the  quarter-chord,  and  the  lowest  peak  has  moved  to  about 
two  chords  behind  the  airfoil.  Also,  a second  peak  has  begun  to  form 
downstream. 

Figure  14  shows  the  velocity  vector  field  at  the  three  times  of  the 
previous  figure.  The  figure  shows  the  change  in  the  angle  of  attack  on 
the  airfoil,  due  to  the  free  surface  movement,  as  time  increases.  The 
velocity  at  the  trailing  edge  shows  that  the  fluid  leaves  the  trailing 

edge  smoothly  creating  a wake  behind  the  airfoil.  The  value  of  conti- 

nuity at  the  trailing  edge  was  the  largest  value  in  the  field.  Due  to 
the  low  Reynolds  number,  the  boundary  layer  is  very  thick. 

Figure  15,  shows  the  pressure  distribution  at  three  time  steps. 

The  figure  shows  the  dominance  of  the  inertial  forces  increasing  over 
the  buoyancy  forces  as  time  increases  from  2 to  6. 

Figures  16  and  17  show  the  effect  of  Reynolds  number  on  the  free 
surface  at  a constant  Froude  number  and  constant  time.  Figure  16  shows 
two  coordinate  systems,  one  at  a Reynolds  number  of  20  and  the  other 
at  a Reynolds  number  of  100.  The  Re  - 100  case  has  less  effect  on  the 
free  surface  than  does  the  Re  = 20  case.  The  effect  of  the  Reynolds 
number  on  the  free  surface  is  due  to  the  smaller  boundary  layer  on  the 

airfoil  for  a Reynolds  number  of  100  than  for  a Reynolds  number  of  20. 

Figure  17  shows  the  pressure  distributions  about  a Karman-Treff tz 
airfoil  for  two  Reynolds  numbers.  Re  = 20  and  Re  = 100.  The  constant 
parameters  are  F » 0.5  and  t = 6.  The  pressure  coefficients  for  Re  = 
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20  are  larger  than  the  pressure  coefficients  for  Re  = 100.  The  lift 
coefficient  for  both  cases  is  due  mainly  to  the  buoyancy  forces.  Figure 
18  shows  a velocity  vector  field  for  a Karman-Tref f tz  airfoil.  Re  = 100, 

F = 0.5,  and  t.  = 6.  The  boundary  layer  is  shown  to  be  smaller  on  the 
airfoil  than  in  the  Re  = 20  case.  The  trailing  edge  wake  is  well 
defined . 

Figures  19  and  20  show  the  versatility  of  the  coordinate  transfor- 
mation. Figure  19  shows  the  coordinate  system  for  a circular  cylinder 
located  one  chord  below  the  free  surface.  Three  times  are  shown.  The 
flow  parameters  are  Re  = 20  and  F = 0.5  at  a time  of  6.  Two  wave  peaks 
are  shown  on  the  free  surface.  The  circular  cylinder  affects  the  free 
surface  more  than  the  airfoil,  which  should  be  expected.  Figure  20  shows 
the  velocity  vectors  for  the  circular  cylinder  of  Figure  19.  At  a time 
of  6,  the  stagnation  point  started  to  move  up  the  front  of  the  cylinder. 

Results  are  presented  in  Figures  21-23  for  a hydrofoil  in  pitching, 
plunging,  and  longitudinal  oscillation.  All  solutions  are  presented  in 
the  free  stream-fixed  coordinate  reference  frame.  Also  all  solutions  were 
run  using  the  flow  parameters.  Re  = 20  and  F = 1.0.  The  same  coordinate 
system  was  used  for  all  solutions.  A Karman-Tref ftz  airfoil  was  placed 
at  one  chord  below  the  free  surface.  The  infinity  boundary  was  located 
10  chords  from  the  airfoil. 

In  Figure  21,  the  airfoil  is  moving  in  the  negative  x-direction.  This 
solution  is  equivalent  to  the  other  solutions  of  previous  sections.  How- 
ever, the  hydrofoil  is  moving  relative  to  the  "outer"  boundary.  The  coor- 
dinate lines  close  to  the  body  are  moving  with  the  airfoil.  The  free  stream- 
fixed  reference  frame  clearly  shows  the  fluid  being  pushed  by  the  airfoil. 
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Also,  fluid  is  moving  in  at  the  trailing  edge  to  fill  the  space  evacuated 
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by  the  airfoil. 

In  Figure  22,  the  airfoil  is  moving  toward  the  free  surface.  As  the 
airfoil  moves  toward  the  free  surface,  fluid  is  pushed  up  and  to  the  sides. 

Vortices  are  created  at  the  leading  edge  and  at  the  trailing  edge  of  the 
airfoil  because  the  fluid  is  moving  from  the  upper  side  to  the  lower  side. 

The  fluid  that  is  moved  up  by  the  airfoil  disturbs  the  free  surfaces  by 
pushing  up  the  free  surface  above  the  airfoil. 

In  Figure  23,  the  airfoil  is  pitching  5°  about  its  center  chord. 

Three  times  (t  = 1,  5,  8)  are  shown.  The  airfoil  takes  a time  of  10  to 
pitch  from  0°  to  5°  and  back  to  0°.  The  vortices  can  be  seen  forming  as 
the  airfoil  pitches. 

Some  of  the  solutions  were  generated  on  the  UNIVAC  1106  single  pro- 
cessor and  the  latest  solutions  were  generated  on  the  upgraded  UNIVAC 
1106  dual  processor.  There  are  many  factors  which  determine  the  com- 
puter time  required  for  a solution,  for  example,  the  way  the  object  pro-  1 

gram  is  loaded  in  the  computer  code.  Several  examples  of  time  required 
to  solve  the  coordinate  system  and  Navier-Stokes  equations  are  pre- 
sented . 

The  uncontracted  coordinate  system  requires  from  3 to  6 minutes  to 
converge  depending  on  the  field  size  and  convergence  criteria.  Depend- 
ing on  the  type  of  attraction  required,  the  contracted  coordinate  sys- 
tem takes  up  to  30  minutes.  j 
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Two  similar  solutions  were  generated  with  field  sizes  of  54  x 30 
and  54  x 60.  The  acceleration  parameter  for  pressure  was  1.8  and  the 
acceleration  parameter  for  velocity  was  .8.  The  constant  of  .1  was 
used  in  the  pressure  iteration  on  the  body.  The  flow  parameters  were 
Re  = 20  and  F = 1.  The  54  x 30  solution  took  239  minutes  to  generate  600 
time  steps  and  the  54  x 60  solution  took  452  minutes  to  generate  600  time 
steps.  The  maximum  number  of  iterations  for  a time  step  to  converge  was 
11  at  time  step  313  for  the  54  x 30  field.  For  the  54  x 60  field  the 
maximum  number  of  iterations  was  12  at  time  step  304. 
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VII.  HYDROFOIL  IN  FREE  SURFACE  SOLUTION  (Ref.  [3]) 


The  physical  and  transformed  planes  used  in  this  solution  are  shown 
in  Figure  24.  Since  the  hydrofoil  is  in  the  free  surface,  rather  than 
submerged,  the  physical  region  is  simply-connected,  its  boundaries  being 
the  wetted  portion  of  the  hydrofoil  contour, (2)  - Q),  the  free  surface 
fore,  O - Q),  and  aft, (3)  - (jj),  and  a remote  semi-circular  boundary, CD  ~ 
0 located  at  a sufficient  distance  from  the  hydrofoil  to  be  undisturbed 
by  the  flow.  The  transformed  plane  is  a rectangle,  with  the  wetted  portion 
of  the  hydrofoil  contour  transforming  to  the  upper  horizontal  side,  the 
free  surface  fore  and  aft  transforming  to  the  left  and  right  vertical  sides, 
respectively,  and  the  semi-circular  remote  boundary  transforming  to  the 
lower  horizontal  side  as  indicated  in  Figure  24.  This  configuration  differs 
from  that  used  for  the  submerged  hydrofoil  in  that  the  physical  region  was 
doubly-connected  with  the  submerged  body. 

As  noted  above,  the  curvilinear  coordinates  (C»h)  are  taken  as  the 
solution  of  two  elliptic  partial  differential  equations.  The  particular 
equations  used  in  the  present  solution  are  those  of  [15],  which  differ  from 
the  original  system  of  [9],  used  for  the  submerged  hydrofoil,  only  in  the 
form  of  the  coordinate  system  control  terms  (the  terms  involving  the  func- 
tions P and  Q below  and  in  Eq . (2)  above.)  Thus  £ and  n are  determined  by 
the  solution  of 

'xx  + «,y  ■ (tx2  - S2)  P<£'n>  <17a) 

"xx  + V ■ <nx2  + ny2>  <!«•">  <1,b) 

With  reference  again  to  Figure  24,  the  boundary  conditions  for  these 
equations  are  as  follows: 
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(a)  on  t*.e  wetted  portion  of  the  hydrofoil  contour,  (5)  - 

n = 12  = constant,  £ varying  monotonically  from  £ to 

S2  U2  > Cx)  from©  to© 

(b)  on  the  free  surface,  (T)  - (©  and©  - ©: 

i = = constant  on©  - © C = C2  = constant  > £ on©  - © 

n varying  monotonically  from  to  n2 
(n2  > 1-^)  from©to©and  from©to© 

(c)  on  the  remote  boundary,©  - (4) 

1=0,=  constant<  n2.  £ varying  monotonically  from 

Cj.  to  i2  from©  to© 

The  coordinate  control  functions,  P(£,n)  and  Q(£,n)  serve  to  concen- 
trate coordinate  lines  as  desired  to  resolve  expected  large  gradients. 

The  theory  and  use  of  such  control  has  been  discussed  in  some  detail  in 
[9].  In  the  present  application,  these  functions  are  determined  from  the 
specified  spacing  of  points  on  the  hydrofoil  contour  and  free  surface, 
those  on  the  body  being  concentrated  near  the  free  surface  and  those  on 
the  surface  being  concentrated  near  the  hydrofoil  as  in  Figure  25.  The 
details  of  this  determination  of  P and  Q are  given  in  Reference  3.  The 
resultant  concentration  of  coordinate  lines  near  the  body  and  free  surface 
is  evident  in  Figure  25. 

In  the  initial  stages  of  this  study,  the  control  functions,  P and 
Q,  were  taken  as  sums  of  decaying  exponentials  that  cause  attraction  of 
coordinate  lines  to  specified  lines  and/or  points  as  used  for  the  submerged 
hydrofoil.  Some  of  the  results  given  below  were  obtained  on  coordinate 
systems  using  this  type  of  control  as  will  be  noted.  The  new  control 
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procedure  has  the  advantage  of  automating  the  control  and  eliminating  the 
need  for  judgmental  estimation  of  the  attraction  amplitudes  and  decay 
factors  necessary  to  achieve  a desired  degree  of  line  concentration. 

With  the  current  modification  in  the  control  functions,  Eq . 2 in 
the  transformed  plane  are  replaced  by 

axK  - 2Bx?n  + Yxnn  + aPx^  + yQx^  = 0 (18a) 

av ££  ~ 28y^n  + yynri  + aPy^  + yQy^  = 0 (18b) 

The  boundary  conditions  for  x and  y are  as  follows: 

(a)  On  the  hydrofoil: 

x and  y specified  by  the  chosen  spacing  of  points  around  the  hydrofoil 
contour.  These  points  move  on  the  contour  with  time  as  a result  of 
motion  of  the  hydrofoil  and  also  because  of  movement  of  the  free 
surface-body  contact  points  on  the  contour. 

(b)  On  the  free  surface: 

Xj-  = 0 initially,  fixed  thereafter,  y from  the  surface  movement, 

Eq.  (10).  (The  first  of  these  allows  the  points  to  slide  along  the 
free  surface  so  that  the  coordinate  lines  are  initially  vertical  at 
the  free  surface.) 

(c)  On  the  remote  boundary: 

x and  y fixed  and  specified  by  the  chosen  spacing  of  points  along 
the  remote  boundary. 

This  point  distribution  on  the  body  and  remote  boundary  was  taken 
according  to  equi-angular  spacing  over  the  body  and  remote  boundary  arcs 
in  the  earlier  stages  of  the  present  investigation.  Later  an  unequal 
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spacing  was  used,  with  the  angular  separation  of  points  varying  on  a 
sine  curve,  so  that  the  closest  spacing  occurs  adjacent  to  the  free  sur- 
face . 

Since  the  free  surface  deforms  in  time,  with  consequent  motion  of 
its  intersections  with  the  hydrofoil  contour,  only  the  relative  dis- 
tribution of  points  on  the  hydrofoil  contour  is  kept  fixed.  The  points 
thus  slide  along  the  wetted  portion  of  the  hydrofoil  while  maintaining 
the  same  relative  spacing  from  adjacent  points  as  time  progresses.  This 
is  accomplished  by  locating  the  points  on  the  contour  at  fixed  percent- 
ages of  the  angle  subtended  by  the  arc  between  the  two  intersections 
with  the  free  surface.  This  subtended  angle,  of  course,  changes  in 
time.  When  the  hydrofoil  oscillates,  the  movement  follows  the 
oscillating  motion  of  the  body  as  well. 

The  points  on  the  free  surface  are  initially  determined  by  a Neumann 
boundary  condition  that  requires  the  coordinate  lines  to  be  vertical  at 
the  moving  surface.  The  local  elevation  of  the  surface  is  determined  by 
the  equations  of  motion  for  the  free  surface  as  discussed  in  Section  III. 

The  points  thus  slide  along  the  free  surface  as  the  surface  deforms  in  time. 

Results  of  the  numerical  solution  are  presented  for  a cir.  ;.lar 
cylinder  hydrofoil  in  two  flow  configurations: 

(a)  Accelerating  translational  motion  parallel  to  the  plane  of 
the  initially  undisturbed  flat  free  surface. 

(b)  Oscillatory  plunging  motion  normal  to  the  plane  of  the 
initially  undisturbed  flat  free  surface. 

In  each  case  the  axis  of  the  cylinder  is  in  the  plane  of  the  initially 

undistrubed  flat  free  surface.  The  fluid  is  physically  unbounded  except 

by  the  free  surface,  with  no  disturbance  remote  from  the  hydrofoil. 


In  the  translational  case  (a)  the  acceleration  is  linear,  with 


the  Reynolds  and  Froude  numbers  given  by 
R = 20 1 


F = 2t 


these  numbers  being  based  on  the  cylinder  diameter  and  current  velocity. 

For  the  plunging  case  the  motion  of  the  hydrofoil  is  sinusoidal 
with  the  elevation  of  the  cylinder  axis  relative  to  the  plane  of  the 
initially  undisturbed  free  surface  given  by 


y = A sin(— ) 

where  A and  P are  the  amplitude  and  period,  respectively,  of  the  motion. 
The  velocity  of  the  cylinder  is  thus 


y 


2ttA 

P 


,2irt 

cosv-^- 


) 


and  the  Reynolds  and  Froude  numbers  are  then  given  by  R = 20y  and  F - 2y , 
respectively . 


In  each  case  the  coordinate  system  is  of  the  form  shown  in  Figure  25 

and  discussed  above  in  this  section,  with  37  points  on  the  body  and 

30  points  on  the  free  surface  on  each  side  of  the  body.  The  remote 

boundary  where  the  fluid  is  undisturbed  is  located  at  a radius  of  10 

cylinder  diameters.  The  convergence  acceleration  parameters  used  were 

1.0  for  the  momentum  equations  (7a-b),  1.8  for  the  Poisson  equation  (7c), 

0.45  for  the  surface  pressure  equation  (9),  and  1.85  for  the  coordinate 

system  equations  (18).  The  iterative  convergence  criteria  used  were  10 

-4 

for  the  coordinate  system  and  10  for  the  velocity  and  pressure. 
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The  initial  point  distribution  on  the  cylinder  was  determined  by 
a sine  curve,  with  points  distributed  symmetrically  and  concentrated  near 
each  free  surface  contact  point.  The  initial  distribution  on  the  fret 
surface  was  determined  by  an  exponential  curve,  with  points  concentrated 
near  the  body.  It  was  found  necessary  to  have  the  points  adjacent  to 
the  free  surface-body  contact  point  approximately  equidistant  from  the 
contact  point  else  stability  problems  arose  with  the  surface.  The  points 
move  on  both  the  hydrofoil  and  free  surface  as  time  passes,  but  the  same 
relative  distributions  are  maintained  as  discussed  in  the  previous  sec- 
tion. 

As  noted  in  the  discussions  above,  the  curvilinear  coordinate  system 

continually  deforms  as  time  progresses,  always  keeping  a coordinate  line 
coincident  with  the  deforming  free  surface.  This  behavior  is  evident  in 
Figure  25  which  shows  the  coordinate  system  at  four  times  for  the  trans- 
lating hydrofoil.  The  free  surface  rises  in  front  of  the  hydrofoil,  and 
the  fore  contact  point  slides  up  along  the  hydrofoil  contour.  At  the  rear 
of  the  hydrofoil,  the  surface  falls,  and  the  aft  contact  point  moves  down- 
ward. 

Velocity  vectors  and  the  hydrofoil  pressure  distribution  for  this 
solution  are  shown  at  one  time  in  Figure  26.  The  vectors  clearly  show 
the  fore  and  aft  stagnation  points  to  be  well  below  the  corresponding  sur- 
face contact  points  on  the  hydrofoil.  The  pressure  distribution  shows 
a positive  pressure  spike  adjacent  to  both  contact  points,  but  . smooth 
distribution  elsewhere  on  the  hydrofoil.  This  spike  is  due  to  numerical 
error  resulting  probably  from  the  modeling  of  the  contact  point  movement. 
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Figures  27  and  28  give  the  surface  elevation  and  pressure  distribution 
in  the  vicinity  of  the  hydrofoil  for  three  times.  There  is  a smooth  rise 
of  the  surface  in  front  of  the  hydrofoil  and  a smooth  depression  behind. 

The  low  Reynolds  number  for  this  solution  causes  distributions  to  dissi- 
pate downstream  so  that  no  surface  waves  occur.  The  surface  pressure 
distribution  in  Figure  28  shows  some  jaggedness  just  upstream  of  the  hydro- 
foil. This  is  probably  due  to  the  numerical  phenomenum  of  wiggles  (two- 
cell wavelength  oscillations),  introduced  in  the  present  case  by  the 
modeling  of  the  contact  point  where  the  pressure  is  fixed  at  zero. 

The  initially  undeformed  coordinate  system  used  in  the  oscillatory 
solution  is  shown  in  Figure  29.  A stronger  concentration  of  lines  near  the 
free  surface  and  hydrofoil  contour  was  used  in  view  of  the  results  dis- 
cussed above.  Figure  30  shows  the  temporal  oscillation  of  the  lift  coef- 
ficient. The  curve  is  seen  to  be  deformed  from  a pure  sinusodial  oscilla- 
tion. After  an  initial  rise,  the  force  remains  upward  throughout  the  cycle. 

Figure  31  shows  a series  of  plots  of  velocity  vectors  at  several 
times  during  the  cycle,  while  Figure  32  shows  the  same  thing  but  in  detail 
of  the  region  around  the  right  surface-body  contact  point.  As  the  body 
rises  initially,  the  fluid  moves  downward  from  the  surface  and  inward 
toward  the  void  being  left  by  the  rising  body  (cf.  T = 0.01  in  Figure  31). 

Since  the  fluid  adjacent  to  the  body  must  rise  with  the  body,  however,  due 
to  the  viscous  no-slip  boundary  condition,  a vortex  is  created  near  the 
contact  point  (T  = O.Or  in  Figure  32).  (A  similar  vortex,  but  of  the 
opposite  rotation  is,  of  course,  created  off  the  left  contact  point.) 

These  vortices  move  away  from  the  body  and  decrease  in  intensity  as  the 
quarter-cycle  is  approached  (cf.  T = 0.11  in  Figure  32). 
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At  the  quarter-cycle  (0.157),  the  body  reaches  its  highest  point  and 
then  reverses  its  motion  to  move  downward.  This  forces  the  fluid  beneath 
the  body  to  the  sides.  The  beginning  of  this  sideward  motion  can  be  seen 
at  T = 0.17  in  Figure  31  just  beneath  the  body.  At  this  time,  inertia 
causes  most  of  the  fluid  to  still  reflect  the  previous  upward  movement  of 
the  body.  This  inertial  effect  is  evident  also  in  the  corresponding 
detail  plot  in  Figure  32,  where  the  fluid  adjacent  to  the  body  has 
reversed  its  motion  and  is  moving  downward  with  the  body  while  the  rest  of 
the  fluid  motion  is  qualitatively  similar  to  that  at  T = 0.15  before  the 
quarter-cycle . 

As  time  passes,  the  influence  of  the  downward  motion  of  the  body 
spreads  progressively  throughout  the  fluid  so  that  more  and  more  of  the 
fluid  acquires  downward  and  sideward  motion  beneath  the  body,  with  con- 
sequent upward  motion  toward  the  surface  (cf.  T = 0.21  and  0.31  in  Figure 
31).  This  annihilates  the  vortices,  and  new  vortices  of  opposite  rotation 
to  the  original  form  just  off  each  surface-body  contact  point  (cf.  T = 

0.21  in  Figure  32.)  These  vortices  also  move  away  from  the  body  and 
decrease  in  intensity  as  the  body  moves  toward  its  lowest  point  at  the 
three-quarter  cycle  time  (T  = 0.471). 

At  this  time  the  motion  of  the  body  again  reverses,  and  the  body 
starts  back  upward.  This  causes  inward  motion  to  begin  just  beneath  the 
body  (T  = 0.49  in  Figure  31)  with  upward  motion  adjacent  to  the  body 
(T  = 0.49  in  Figure  32).  This  new  pattern  of  motion  then  spreads  out  into 
the  remainder  of  the  fluid,  competing  initially  with  the  inertial lv  per- 
sisting motion  from  before  the  last  body  reversal.  As  at  the  quarter- 
cycle,  the  existing  vortices  are  annihilated,  and  new  ones  of  opposite 
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rotation  again  form  off  the  contact  points  (cf.  T = 0.53  in  Figure  32). 
The  general  fluid  motion  is  again  downward  from  the  surface,  with  inward 
and  upward  motion  beneath  the  rising  body  (cf.  T = 0.61  in  Figure  31)  as 


at  the  beginning  of  the  cycle. 

The  movement  of  the  hydrofoil  free  surface  contact  points  along  the 
hydrofoil  contour  is  modeled  by  a condition  of  continuity  as  discussed  in 
detail  in  Ref.  3.  Essentially  this  model  causes  the  contact  points  to  slide 
along  the  hydrofoil  contour  in  response  to  a net  inbalance  of  flow  into  the 
cell  at  the  contact  point.  Net  inflow  will  thus  cause  that  contact  point 
to  slide  upward  along  the  contour.  Another  model  based  on  a condition  of 
zero  stress  at  the  contact  point  was  also  investigated  but  was  found  to  be 
unsatisfactory  as  also  discussed  in  Ref.  3. 

Figures  33  and  34  show  the  surface  elevation  and  pressure  at  approxi- 
mately the  quarter,  half,  three-quarter , and  full  cycle  times.  The  ele- 
vation curves  show  that  the  mean  surface  position  is  not  flat,  but  is 
depressed  in  the  vicinity  of  the  body.  This  result  is  in  qualitative  agree- 
ment with  a periodic  boundary  layer  solution  and  experimental  flow  visuali- 
zation results  given  in  Schlichting  [14]  for  a circular  cylinder  oscillating 
in  an  unbounded  fluid.  There  it  is  shown  that  a mean  secondary  motion  exists 
in  which  fluid  moves  from  the  sides  toward  the  body  (normal  to  the  direction 
of  oscillation)  and  then  away  from  the  body  parallel  to  the  oscillation 
direction  (cf.  Figure  11.7  of  [14]).  In  the  present  case  this  type  of  mean 
flow  would  be  toward  the  body,  parallel  to  the  free  surface,  and  then  away 
from  the  surface  beneath  the  body.  This  then  would  result  in  a mean  surface 
depression. 


38 


VIII.  CONCLUSION 


The  excellent  comparison  of  the  numerical  results  based  on  the  techni- 
que of  numerically-generated  boundary- fitted  coordinate  systems  with  the 
Blasius  boundary  layer  solution  for  the  semi- infinite  flat  plate  demon- 
strates that  this  technique  can  be  quite  useful  in  the  numerical  solution 
of  the  Navier-Stokes  equations. 

The  technique  of  numerically  generated  boundary-fitted  coordinate 
systems  is  clearly  an  effective  aid  in  treating  flow  problems  involving 
both  free  surfaces  and  solid  boundaries.  With  this  technique  the  complica- 
tion of  the  boundary  shape  is  essentially  removed  from  the  problem.  It  is 
possible  to  obtain  numerical  solutions  for  viscous  flow,  with  viscous 
boundary  conditions  on  the  free  surface  as  well  as  on  the  solid  body. 

The  research  results  presented  in  this  report  leave  several  problems 
unresolved.  Regarding  the  submerged  hydrofoil  solution,  the  coordinate 
configuration  used  had  a zero  Jacobian  between  grid  points  in  the  field  of 
calculation.  This  zero  Jacobian  made  it  difficult  to  contract  coordinate 
lines  near  the  branch.  Also,  an  ambiguity  in  the  finite  difference  expres- 
sions for  the  cross  derivatives  at  the  point  of  zero  Jacobian  led  to 
ambiguous  results  in  the  coordinate  solutions.  Thirdly,  the  coordinate  con- 
trol functions  were  found  to  be  inadequate  in  controlling  coordinate  lines 
in  the  field.  An  arbitrary  change  in  the  coordinate  control  often  led  to 
unpredictable  results  for  the  physical  coordinates.  Some  progress  was  made 
in  this  area  during  the  latter  stages  concerning  the  hydrofoil  in  the  free 
surface  with  the  incorporation  of  an  automated  control.  Finally,  at  pro- 
ject termination,  the  solution  could  not  be  run  successfully  for 
Reynolds  numbers  greater  than  100  for  the  submerged  hydrofoil.  The  combina- 
tion of  the  coordinate  contraction  functions  and  the  field  point  witli  a zero 
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Jacobian  is  believed  to  be  the  cause  of  a pressure  source  that  occurred  at 
the  trailing  edge  of  hydrofoils  for  Re  > 100.  Finally,  since  higher 
Reynolds  numbers  were  not  obtained,  the  method  could  not  be  verified  with 
the  experimental  data  available  for  submerged  hydrofoils. 

The  presence  of  a zero  Jacobian  in  the  field  is  not  a universal  fea- 
ture in  the  boundary-fitted  coordinate  systems,  but  is  peculiar  to  the 
type  of  configuration  adopted  for  the  transformed  plane  in  the  submerged 
hydrofoil  solution.  This  configuration  did  have  certain  advantages  as 
noted  in  spite  of  the  presence  of  the  zero  Jacobian.  The  configurations 
used  for  the  hydrofoil  in  the  free  surface  and  for  the  semi-infinite  flat 
plate  (and  for  the  external  flow  about  airfoils  in  other  studies)  do  not 
have  any  zeros  of  the  Jacobian  in  the  field.  Further  study  would  be 
necessary  to  develop  better  configurations  for  the  submerged  hydrofoil  case. 

Concerning  the  hydrofoil  in  the  surface,  the  coordinate  configuration 
was  less  of  a problem,  and  no  zeros  of  the  Jacobian  occurred  in  the  field. 

The  results  given  in  the  present  work  are  all  at  very  low  Reynolds  number, 
but  the  solution  can  in  principle  be  run  at  any  Reynolds  number  by  increasing 
the  attraction  of  the  coordinate  lines  to  the  body  and  free  surface  at 
higher  Reynolds  numbers  in  order  to  maintain  a sufficient  number  of  lines 
in  the  viscous  layers.  Such  a procedure  is  currently  under  investigation 
in  connection  with  the  flow  about  airfoils.  The  problem  is  made  more 
difficult,  however,  with  increasing  Reynolds  number.  More  investigation 
of  the  control  of  the  coordinate  system  so  that  sufficiently  close  spacing 
is  maintained  near  the  free  surface  as  it  deforms  is  necessary,  as  is 
further  study  of  the  modeling  of  the  hvdrofoil-free  surface  contact  point 
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movement . 


The  inclusion  of  a bottom  of  arbitrary  shape  in  the  current  solution 


is  not  difficult,  nor  would  be  the  addition  of  a second  hydrofoil  in  tandem. 
(The  inclusion  of  a solid  bottom  was  accomplished  in  an  extension  of  the 
present  work  as  reported  in  [16].)  Similarly  the  hydrofoil  and/or  bottom 
could  be  allowed  to  deform  in  time  without  complicating  the  problem 
unduly.  This  is  because  all  of  the  computation  is  done  on  the  fixed 
rectangular  transformed  grid  regardless  of  the  shape  or  movement  of  the 
physical  boundaries.  Wind  shear  on  the  free  surface  could  also  be  added 
by  a change  in  surface  boundary  conditions  to  include  applied  external 
shear  as  well  as  pressure. 


41 


REFERENCES 


1.  Walker,  R.  L.,  "Numerical  Solution  of  the  Navier-Stokes  Equations 
for  Incompressible  Viscous  Laminar  Flow  Past  a Semi-Infinite  Flat 
Plate,”  M.S.  Thesis,  Mississippi  State  University,  Mississippi 
State,  Mississippi  (1974). 

2.  Shanks,  S.  P.,  "Numerical  Simulation  of  Viscous  Flow  about  Submerged 
Arbitrary  Hydrofoils  using  Non-orthogonal , Curvilinear  Coordinates," 
Ph.D.  Dissertation,  Mississippi  State  University,  Mississippi  State, 
Mississippi  (1977). 

3.  Thompson,  J.  F.,  and  Shanks,  S.  P.,  "Numerical  Solution  of  the 
Navier-Stokes  Equations  for  2D  Surface  Hydrofoils,"  EIRS-ASE-77-4 , 
Mississippi  State  University,  Mississippi  State,  Mississippi,  (1977). 

4.  Thames,  F.  C.,  Thompson,  J.  F. , et . al . , "Numerical  Solutions  for 
Viscous  and  Potential  Flow  about  Arbitrary  Two-Dimensional  Bodies 
using  Body-Fitted  Coordinate  Systems,"  accepted  for  publication  in 
Journal  of  Computational  Physics  (1976). 

5.  Thompson,  J.  F. , Thames,  F.  C.,  Walker,  R.  L. , Shanks,  S.  P., 
"Numerical  Solutions  of  the  Unsteady  Navier-Stokes  Equations  for 
Arbitrary  Bodies  using  Boundary-Fitted  Curvilinear  Coordinates," 
Proceedings  of  Arizona/AFOSR  Symposium  on  Unsteady  Aerodynamics, 
University  of  Arizona  (1975). 

6.  Thompson,  J.  F.,  Thames,  F.  C.,  Shanks,  S.  P.,  Hodge,  J.  K. , and 
Mastin,  C.  W. , "Solutions  of  the  Navier-Stokes  Equations  in  Various 
Flow  Regimes  on  Fields  Containing  any  Number  of  Arbitrarily-Shaped 
Two-Dimensional  Bodies  using  Boundary-Fitted  Curvilinear  Coordinate 
Systems,"  5th  International  Conference  on  Numerical  Methods  in 
Fluid  Dynamics,  Enschede,  the  Netherlands  (1976). 

7.  Thompson,  J.  F. , Thames,  F.  C.,  and  Mastin,  C.  W. , "Automatic 
Numerical  Generation  of  Body-Fitted  Curvilinear  Coordinate  System 
for  Field  Containing  any  Number  of  Arbitrary  Two-Dimensional 
Bodies,"  Journal  of  Computational  Physics,  15,  299  (1974). 

8.  Thompson,  J.  F. , Thames,  F.  C.,  Mastin,  C.  W. , "TOMCAT  - A Code  for 
Numerical  Generation  of  Boundary-Fitted  Curvilinear  Coordinate 
Systems  on  Fields  Containing  any  Number  of  Arbitrary  Two-Dimensional 
Bodies,"  accepted  for  publication  in  Journal  of  Computational 
Physics  (1976). 

9.  Thompson,  J.  F.,  Thames,  F.  C.,  and  Mastin,  C.  W. , "Boundary- 
Fitted  Coordinate  Systems  for  Solution  of  Partial  Differential 
Equations  on  Fields  Containing  any  Number  of  Arbitrary  Two- 
Dimensional  Bodies,"  NACA  CR-2729,  1977. 


42 


10.  Hodge,  J.  L.,  "Numerical  Solution  of  Incompressible  Laminar  Flow 
about  Arbitrary  Bodies  in  Body-Fitted  Curvilinear  Coordinates," 

Ph.D.  Dissertation,  Mississippi  State  University,  Mississippi  State, 
Mississippi  (1975). 

11.  Hirt,  C.  W.,  and  Harlow,  F.  H. , "A  General  Corrective  Procedure  for 
the  Numerical  Solution  of  Initial  Value  Problems,"  Journal  of 
Computational  Physics,  2,  114  (1967). 

12.  Thames,  F.  C.,  "Numerical  Solution  of  the  Incompressible  Navier- 
Stok.es  Equations  about  Arbitrary  Two-Dimensional  Bodies,"  Ph.D. 
Dissertation,  Mississippi  State  University,  Mississippi  State, 
Mississippi  (1975). 

13.  M.  Israeli,  "A  Fast  Implicit  Method  for  Time  Dependent  Viscous 
Flows,"  Studies  in  Applied  Math,  49,  327  (1970). 

14.  Schlichting,  H. , Boundary  Layer  Theory,  4th  Ed.,  McGraw-Hill,  (1960). 

15.  Warsi,  Z.  U.  A.,  Thompson,  J.  F. , "Machine  Solutions  of  Partial 
Differential  Equations  in  the  Numerically  Generated  Coordinate 
System,"  MMSU-EIRS-ASE-77-1,  Engineering  and  Industrial  Research 
Station,  Mississippi  State  University ,( 1 977) . 

16.  Thompson,  J.  F.  and  Shanks,  S.  P.,  "Numerical  Solution  of  the 
Navier-Stokes  Equations  for  Arbitrary  21)  Submerged  Hydrofoils 
in  a Channel  of  Finite  Depth."  MMSU-EIRS-ASE-77-6 , Mississippi 
State  University,  Mississippi  State,  Mississippi  (1977). 


43 


Transformed  Plane 


Figure  1.  Field  Transformation  - Single  Body 


Physical  Field 


n 


Figure  2 


Transformed  Field 


Relation  between  Physical  and  Transformed 
Fields  - Semi-Infinite  Flat  Plate 


(b) . Actual  System 


. _ 


(a),  x = 0.001  (b).  x - 0.01 


(c).  x = 0.50  (d).  x *=  1.00 


Figure  4.  Semi-Infinite  Flat  Plate  Solution.  t = 1.24  (Downstream 

Boundary  at  Eight  Times  the  Distance  where  Blasius  Boundary 
Layer  Slope  is  0.01,  XMAX  «»  0.45312) 
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Pressure  distribution  for  various  Reynolds  and 
Froude  numbers  at  t = 8.0  - Karman-Tref f tz 
Airfoil  located  1 chord  below  free  surface. 
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Figure  19.  Coordinate  system  at  three  times, 
circular  cylinder  located  1 chord 
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Figure  22.  Velocity  vector  field  of  Karman-Treff tz  Airfoil  translating  toward  free  surface,  Re 

F = 1.0  and  t = 1.75  (free  stream-fixed  coordinates). 


J. 


Velocity  vector  field  of  Karman-Tref f tz  slowly 
pitching  5°  about  its  center  chord  at  three 
times.  Re  = 20,  F = 1.0  (free  stream- fixed 
coordinates) . 
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25.  Deforming  Coordinate  System  - Translating  Hvdrofoil 
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Figure  27.  Surface  Elevation  - Translating  Hydrofoil 
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Figure  31.  Velocity  Vectors  - Plunging  Hydrofoil  (0.001  Amplitude,  0.628  Period) 
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